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The  spherical  agglomerate  model  represents  the  most  detailed  description  of  the  PEM  fuel  cell  catalyst 
layer  as  it  accounts  for  both  micro  and  macroscale  transport  phenomena.  The  usual  approach  with  the 
classical  spherical  agglomerate  model  is  to  couple  the  homogenous  mixture  assumption  for  the 
agglomerate  core  to  its  idealized  spherical  geometry  to  obtain  an  analytical  solution  which  is  easily 
incorporated  within  a  macroscale  model.  In  this  study,  we  incorporate  numerical  results  from  a  modified 
agglomerate  model  based  on  discrete  platinum  particles  [33]  to  create  a  more  physically  realistic  3D 
macroscale  model  for  the  PEM  fuel  cell  cathode  catalyst  layer.  Results  from  the  3D  cathode  model  based 
on  the  modified  particle  approach  are  compared  with  the  classical  model  and  the  macro-homogenous 
model.  We  find  that,  similar  to  the  classical  approach,  the  modified  3D  model  is  able  to  reproduce 
previously  reported  trends  for  reactant,  reaction  rate,  and  overpotential  distributions,  whereas  the 
macro-homogenous  model  fails  to  predict  mass  transport  limitations  properly.  It  is  also  shown  that, 
unlike  the  classical  approach,  the  modified  3D  model  is  able  to  accurately  predict  the  effect  of  Pt  loading 
in  the  diffusion-loss  region. 

©  2013  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Polymer  electrolyte  membrane  (PEM)  fuel  cells  are  one  of  the 
most  environmentally-friendly  energy  conversion  devices,  espe¬ 
cially  for  the  automotive  industry.  During  the  last  decade,  experi¬ 
mental  and  numerical  approaches  have  been  employed  to 
investigate  various  aspects  of  PEM  fuel  cells.  As  a  part  of  that 
ongoing  research,  computational  models  play  an  important  role  in 
the  design  and  optimization  of  PEM  fuel  cells.  The  fidelity  of 
computational  models  is  highly  dependent  on  an  accurate 
description  of  the  electrochemical  phenomena  occurring  in  the 
catalyst  layer  (CL)  due  to  their  critical  impact  on  the  performance 
and  efficiency  of  the  fuel  cell. 
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It  is  a  challenge  to  accurately  model  the  CL  because  its  structural 
and  functional  properties  span  the  range  from  microscopic  to 
macroscopic  scales.  The  complex  porous  CL  structure  containing  Pt, 
carbon  support,  ionomer  electrolyte,  gas  pores,  and  liquid  water 
clearly  represents  a  multiscale,  multiphase  problem  hosting  a  va¬ 
riety  of  reactions,  and  mass  and  heat  transport  phenomena.  For 
device-level  macroscale  simulations,  the  literature  presents  some 
common  approaches  which  can  be  broadly  classified  as  interface 
models  [1,2],  homogeneous  models  [2—5],  thin  film  models  6—9], 
and  agglomerate  models  [10-33]. 

The  simplest  description  of  the  CL  is  provided  by  the  interface 
model  in  which  a  CL  of  infinitesimal  thickness  is  considered  be¬ 
tween  the  membrane  and  the  gas  diffusion  layer  (GDL).  The 
structure  of  the  CL  is  ignored.  The  influence  of  catalytic  activity  is 
introduced  by  source  terms  that  are  used  as  the  boundary  condition 
at  the  membrane/GDL  interface.  This  simplification  neglects  the 
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ohmic  and  diffusion  resistances  of  the  CL.  Therefore,  although  the 
interface  model  is  computationally  efficient,  it  cannot  predict  the 
cell  performance  correctly.  The  3D  PEM  fuel  cell  model  of  Berning 
et  al.  1  ]  represents  a  good  example  of  the  interface  model.  Harvey 
et  al.  [2]  reported  that  current  densities  are  overestimated  by  the 
interface  model. 

The  macro-homogeneous  model  represents  the  CL  as  a  porous 
matrix  consisting  of  a  uniform  mixture  of  Pt,  carbon,  and  electro¬ 
lyte.  Gas  phase  reactants  diffuse  through  the  porous  matrix  to  the 
reaction  sites.  Ions  and  electrons  are  assumed  to  travel  through  the 
electrolyte  and  the  carbon  network,  respectively.  Therefore,  diffu¬ 
sive  and  ohmic  transport  resistances  are  taken  into  account  based 
on  effective  medium  theory.  However,  the  macro-homogeneous 
model  cannot  predict  mass  transport  correctly,  and  is  therefore 
unable  to  provide  a  limiting  current  density.  This  shortcoming  can 
be  clearly  observed  in  the  results  of  Rho  et  al.  [3],  Wang  et  al.  [4], 
and  Eikerling  et  al.  [5]. 

The  thin  film  model  is  almost  identical  to  the  macro- 
homogeneous  model  and  employs  the  same  type  of  CL  structure. 
However,  in  contrast  to  the  macro-homogeneous  model,  the  gas 
pores  are  assumed  to  be  flooded  with  liquid  water.  Therefore, 
reactant  gases  must  first  dissolve  in  liquid  water  in  order  to  be 
transported  to  the  reaction  sites.  This  model  results  in  high  diffu¬ 
sive  resistances  so  a  thin  CL  (i.e.  thin  film)  must  be  assumed  to 
obtain  a  match  with  the  experimental  results.  Unlike  the  homo¬ 
geneous  model,  the  thin  film  model  is  able  to  predict  the  limiting 
current  density  [6-9]. 

Among  the  device-level  modeling  approaches,  the  most  detailed 
representation  of  the  CL  is  provided  by  the  agglomerate  models. 
Unlike  the  previously  mentioned  methods,  only  the  agglomerate 
approach  considers  the  CL  as  a  multi-scale  problem.  The  porous 
structure  of  the  CL  is  assumed  to  consist  of  gas  pores  and  agglom¬ 
erates  (with  an  assumed  microstructure)  that  are  covered  by  an 
electrolyte  film.  The  agglomerate  itself  is  assumed  to  consist  of  a 
homogeneous  mixture  of  carbon-supported  catalyst  and  electrolyte 
as  depicted  in  Fig.  1.  In  this  model,  reactant  gases  travel  through  the 
gas  pores  until  they  reach  the  electrolyte  film  that  covers  the  ag¬ 
glomerates.  The  gases  then  dissolve  into  the  electrolyte  film  and  are 
transported  by  diffusion  inside  the  agglomerate  to  the  reaction  sites. 

In  recent  years,  there  have  been  many  developments  in  CL 
modeling.  Inexpensive  computational  power  has  made  it  feasible  to 
model  the  actual  CL  microstructure  [34-43].  In  addition,  some 
researchers  have  proposed  modifications  to  the  classical  agglom¬ 
erate  model.  For  example,  Jain  et  al.  [30]  compared  spherical,  cy¬ 
lindrical,  and  planar  agglomerate  geometries  and  showed  that  the 
nature  of  the  microscale  geometry  affects  O2  consumption.  Epting 


and  Litster  [31  ]  examined  the  importance  of  agglomerate  size  dis¬ 
tribution  based  on  nano-CT  imaging  of  a  typical  CL,  and  concluded 
that  the  size  distribution  has  a  significant  impact  on  the  results. 
Kamarajugadda  and  Mazumder  [32]  considered  two  overlapped 
agglomerates  with  unequal  radii  as  a  generalized  shape  for 
agglomeration.  They  incorporated  numerical  simulations  of  their 
3D  agglomerate  as  a  sub-grid-scale  model  into  their  2D  PEM  fuel 
cell  model,  and  found  that  while  the  agglomerate  shape  is  unim¬ 
portant  for  small  agglomerates,  it  becomes  significant  for  larger 
agglomerates. 

One  of  the  main  assumptions  of  the  classical  agglomerate  model 
is  that  the  agglomerate  core  consists  of  homogenous  mixture  of 
C|Pt  particles  and  the  ionomer.  A  consequence  of  this  assumption  is 
that  the  model’s  results  become  insensitive  to  Pt  loading  at  high 
current  densities.  Cetinbas  et  al.  [33]  proposed  a  modified  micro¬ 
structure  for  the  agglomerate  core  by  introducing  discrete  Pt  par¬ 
ticles.  They  concluded  that  the  presence  of  discrete  catalyst 
particles  leads  to  a  proper  accounting  of  the  diffusion  process  inside 
the  agglomerate  core.  Therefore,  unlike  the  classical  approach,  the 
modified  model  was  shown  to  accurately  predict  the  effect  of 
catalyst  loading  in  the  diffusion-loss  region. 

In  this  study,  the  microscale  approach  proposed  in  Ref.  33]  is 
incorporated  into  a  3D  single  channel  PEM  fuel  cell  cathode  model. 
However,  unlike  the  classical  agglomerate  approach,  the  modified 
approach  does  not  lend  itself  to  an  analytical  solution  which  is 
easily  coupled  to  the  macroscale  model.  Therefore,  the  coupling 
method  used  in  Ref.  19]  has  to  be  modified  in  order  to  interface 
with  the  numerical  results  from  the  discrete  particle  approach. 
Accordingly,  the  numerical  solution  obtained  from  the  2D  micro¬ 
scale  model  is  expressed  as  a  function  of  macroscale  model  vari¬ 
ables  and  then  incorporated  within  the  3D  cathode  model.  In 
addition,  we  compare  the  predictions  of  the  macro-homogenous, 
classical,  and  modified  agglomerate  approaches  by  plotting  their 
polarization  curves  and  current  density  distributions.  Finally,  we 
investigate  the  sensitivity  of  the  different  agglomerate  models  to  Pt 
loading. 

2.  Cathode  model 

We  have  implemented  a  3D,  single  channel,  isothermal,  steady- 
state  cathode  model  that  accounts  for  the  conservation  of  mass, 
momentum,  chemical  species,  electrons,  and  ions.  The  CL  is 
modeled  by  the  macro-homogenous,  classical,  and  modified 
agglomerate  approaches  which  will  be  explained  in  the  subsequent 
sections.  Symmetry  in  the  single  channel  is  exploited  to  reduce  the 
computational  domain  such  that  it  encompasses  half  of  the  channel 


Ionomer  film 


Fig.  1-  Schematic  of  the  catalyst  layer  structure  in  the  classical  spherical  agglomerate  approach. 
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Fig.  2.  Computational  domain  for  the  3D  cathode  model. 


and  half  of  the  land  area  as  illustrated  in  Fig.  2.  The  domain  di¬ 
mensions  are  listed  in  Table  1.  The  highly-coupled  equations  that 
govern  the  electrochemical  and  transport  phenomena  within  the 
cell  are  solved  with  the  finite  element  software  COMSOL  4.2. 

Assuming  laminar,  incompressible  flow,  mass  and  momentum 
conservation  in  the  gas  channel  are  governed  by  the  following 
equations: 

V  u  =  0  (1) 


Eq.  (5)  solves  for  the  mass  fraction  (w/)  of  each  chemical  species 
with  Rj  representing  the  source  term  due  to  electrochemical  re¬ 
actions  in  the  CL.  The  diffusion  coefficients  for  different  operating 
conditions  can  be  calculated  based  on  the  reference  values  of  binary 
diffusivities  [44]: 

D*-D"'(Pf)(T&)'75  (6) 


p(u-v)u  =  -Vp  +  V-p(Vu  +  (Vu)r)  (2) 

where  u  is  the  velocity  field,  p  is  pressure,  p  is  density  and  p  is 
dynamic  viscosity.  In  the  porous  domains  (GDL  and  CL)  the  flow  is 
described  by  the  Brinkman  equation  which  extends  Darcy’s  law  to 
account  for  the  momentum  balance. 


PV-U  =  Q  (3) 

^((u-v)S|  =  -vp  +  v-^(vu  +  (vu)T)  -  v-^(v-u) 

-(£+«)»  <4> 


where  the  reference  pressure  (pref)  is  1  atm,  the  reference  tem¬ 
perature  (Tref)  is  298  K,  and  the  reference  diffusivities  D^f  used  are 
listed  in  Table  2. 

The  density  of  the  mixture,  which  is  also  used  in  the  mass  and 
momentum  conservation  equations,  is  calculated  using  the  ideal 
gas  law  as  follows. 


P 


MnP 

RT 


where  Mn 


(7) 


The  Bruggeman  approximation  is  commonly  used  to  calculate 
the  effective  properties  of  isotropic  porous  structures.  Here,  we  use 
this  approximation  to  calculate  the  effective  diffusivities. 


where  Q.  is  the  mass  source  term  which  is  calculated  from  the  re¬ 
actions  in  the  CL,  e  is  the  porosity,  and  k  is  the  permeability  of  the 
GDL. 

The  diffusion  of  oxygen,  water  vapor,  and  nitrogen  in  the  cath¬ 
ode  is  described  by  using  the  Maxwell— Stefan  equations  including 
the  effect  of  convection. 


Djf  =  Djke}-5  (8) 

where  q  represents  the  porosity  of  the  corresponding  domain. 
Charge  transport  by  the  electron  and  ion  phases  is  governed  by  the 
following  equations. 

V-(<ffV0m)  =  Vi  (9) 


V-(ffffV0s)  =  -Vi 


(10) 


Table  1 

Geometric  parameters  for  the  single  channel  computational  domain. 


Channel  length 

tch 

50  [mm] 

Table  2 

Channel  height 

bch 

1  [mm] 

Binary  diffusivities  at  1  atm  and  298  I<  [441. 

Channel  width 

Wch 

1  [mm] 

Land  width 

W\ 

1  [mm] 

Do2n2 

2.07  x 

GDL  thickness 

£gdl 

250  [pm] 

doLh2o 

2.64  x 

CL  thickness 

tCL 

15  [pm] 

^n2h2o 

2.64  x 
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where  a^f  and  of ff  are  the  effective  ion  and  electron  conductivities, 
0m  and  0S  are  ionic  and  electronic  potentials,  and  V.i  is  the  volu¬ 
metric  current  density.  The  Bruggeman  approximation  is  again 
used  to  calculate  the  effective  conductivities. 

°f  =  Ojef  (11) 

The  ionic  conductivity  Om)  is  dependent  on  the  local  water 
content  and  temperature.  In  order  to  determine  the  local  ionic 
conductivity  we  used  the  empirical  relation  reported  by  Springer 
et  al.  [45]  for  Nation™. 


(Tm  =  100(5.139  X  10“3A-  3.26  X  10“3)exp 

(12) 

where  X  is  the  local  water  content  which  is  expressed  as  a  function 
of  local  relative  humidity  as  [19]: 


1268 


1 

303 


X  —  0.3  +  1O.80RH  —  16  0rH  +  14.1 


(13) 


where  0 rh  is  the  local  relative  humidity. 

The  effects  of  the  electrochemical  reaction  on  charge  and  reac¬ 
tant  transport  are  described  by  the  source  terms  in  the  governing 
equations.  The  O2  and  H2O  source  terms  in  the  CL  domain  ( R0l  and 
Rh2o)  are  calculated  from  Faraday’s  law: 


R  Vl' 

Ro2  = 


(14) 


rh2o  = 


(^OD  +  1)V-I 
2F 


(15) 


where  «od  is  the  osmotic  drag  coefficient  which  can  be  written  as  a 
function  of  the  nominal  cathode  overpotential  (??nco)  as  [19]: 


1.0 


^nco  <  0-25 


-  ,/1NLU  \ 

^od  —  {  ^Gj^nco  —  31.52t7NCo  +  5.7  0.25  <  ?7nco  —  0 -35 
[  0.3  Vnco  >  0.35 


2  A.  Spherical  agglomerate  model 

Catalyst  layer  modeling  techniques  are  required  to  obtain  an 
expression  for  the  volumetric  current  density  (V.i)  which  quantifies 
the  rate  of  the  electrochemical  reaction  occurring  within  the  CL.  As 
stated  earlier,  the  spherical  agglomerate  approach  is  the  most 
comprehensive  description  of  the  CL  as  it  accounts  for  reaction- 
diffusion  phenomena  inside  the  imaginary  microstructures  called 
agglomerates.  It  is  a  multi-scale  approach  that  ranges  from  the 
agglomerate  at  the  microscale  to  the  entire  CL  at  the  macroscale. 

The  mass  balance  of  O2  in  the  CL  can  be  written  in  terms  of  the 
agglomerate  flux  or  total  reaction  rate  in  the  CL  as  [19]: 

V-No2+aaggN£f  =  0  (20) 

V-iV02  +  =  0  (21) 

where  aagg  is  the  ratio  of  total  agglomerate  surface  area  to  CL  vol¬ 
ume,  Na0f  is  the  O2  flux  at  the  surface  of  a  single  agglomerate,  and 
f?Q2tal  is  the  total  reaction  rate  in  CL.  Sun  et  al.  [19]  tied  the  micro¬ 
scale  agglomerate  to  the  macroscale  CL  by  noting  the  equality  be¬ 
tween  Eqs.  (20)  and  (21)  and  obtained  an  expression  for  V.i. 
Kamarajugadda  and  Mazumder  [32]  used  only  Eq.  (21)  to  integrate 
their  3D  numerical  agglomerate  into  a  2D  device-scale  simulation. 
In  this  study,  we  use  Eq.  (20)  to  derive  the  volumetric  current 
density.  The  goal  is  to  develop  a  procedure  to  couple  the  2D  nu¬ 
merical  results  of  the  modified  agglomerate  problem  in  Ref.  [33] 
with  a  3D  cathode  model. 

The  classical  agglomerate  problem  can  be  summarized  as  a 
reaction-diffusion  process  in  the  homogenous  agglomerate  core 
following  the  pure  diffusion  of  reactant  through  the  ionomer  film. 
For  simplicity,  the  problem  is  solved  in  a  1 D  spherically-symmetric 
domain  as  in  Fig.  3.  The  goal  is  to  calculate  the  reaction  rate  per  unit 
volume  by  obtaining  the  analytical  solution  of  Eqs.  (22)  and  (23). 

(r2^)  =  ^  in  agglomerate  core  (22) 

^37  (r2~ar~)  in  ionomer  film  (23) 


The  nominal  cathode  overpotential  (tjnco)  can  be  defined  as  the 
measure  of  total  losses  (including  activation  and  ohmic  losses)  in 
the  cathode.  It  can  be  expressed  as  19]: 


^NCO  —  0s  |  Current  Collector  0m  I  Membrane- CL  Interface 


(17) 


where  D  is  the  diffusivity  of  O2  in  the  ionomer  film  which  is  curve- 
fitted  in  Ref.  [19]  using  the  experimental  data  of  Ref.  [47]  as: 


D  =  0.0438  exp 


-25  kj  mol 

rT 


-1 


?  -1 

CITT  S  1 


(24) 


where  the  electron  potential  at  the  current  collector  and  ion 
potential  at  the  membrane-CL  interface  are  the  boundary  con¬ 
ditions  for  the  charge  conservation  equations.  All  the  boundary 
conditions  are  clarified  later.  For  the  sake  of  simplicity  our  model 
considers  only  the  cathode  electrode  as  the  majority  of  the 
activation  loss  occurs  in  the  cathode,  whereas  the  reactions  at 
anode  are  several  orders  of  magnitude  faster.  Neglecting  the 
losses  in  anode  and  membrane,  the  cathode  potential  can  be 
expressed  as  [19]: 


Deff  is  the  effective  diffusivity  of  O2  in  the  agglomerate  core 
which  is  calculated  from  the  Bruggeman  correlation,  and  kc  is  the  is 
the  reaction  rate  constant  given  by  the  Butler-Volmer  type 
expression  as: 


kc 


apt 

4F  Q)2,ref 


exp 


exp 


Al  -  ac)F 

v  RT 


(25) 


Vc  =  Eth  -  ^NCO 


Eth  =  1.229-8.456  x  10“4(r- 295.15) +4.31 


x  icr5r 


MPhJ  A+PoJ 


(18) 


(19) 


» 


Agglomerate  Core 


=  1 

1 


Film 


cr 


=  ( 

1 


r  =  0 


i  +  8 


Fig.  3.  ID  spherical  agglomerate  domain  and  boundary  conditions. 
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where  apt  is  the  total  reaction  area  per  unit  volume  of  agglomerate, 
igef  is  the  reference  exchange  current  density,  C02  ref  is  the  standard 
reference  O2  concentration,  ac  is  the  charge  transfer  coefficient,  and 
cp  is  the  local  over-potential  in  the  CL  domain.  The  total  reaction 
area  per  unit  volume  of  the  agglomerate  is  expressed  as: 


apt  =  A) 


mpt 


fa(l  -  ^cl) 


(26) 


where  A0  is  the  catalyst  surface  area  per  unit  mass,  mpt  is  the 
catalyst  loading  (kg  m-2),  and  ta  is  the  thickness  of  the  catalyst 
layer. 

Parthasarathy  et  al.  [46—48]  reported  two  Tafel  slopes  corre¬ 
sponding  to  high  and  low  current  regimes.  In  order  to  account  for 
the  variation  in  Tafel  slope,  we  used  two  distinct  values  for  the 
exchange  current  density  (igef),  and  the  charge  transfer  coefficient 
(ac)  as  given  in  Table  3. 

It  is  important  to  note  that  the  overpotential  within  the 
agglomerate  (<pagg)  is  assumed  to  be  constant  due  to  the  agglom¬ 
erate’s  small  size.  On  the  other  hand,  the  macroscale  overpotential 
(<p)  is  a  variable  governed  by  the  charge  transport  equations. 
Selecting  the  oxygen  electrode  as  the  reference  electrode  with  zero 
potential,  the  local  overpotential  in  the  CL  domain  can  be  expressed 
simply  as  the  difference  between  the  electronic  and  ionic 
potentials. 


(p  —  0s  —  0m 


(27) 


where  Th  is  the  dimensionless  Thiele  modulus  defined  as  the  ratio 
of  the  surface  reaction  rate  and  the  diffusion  rate  through  the 
catalyst  pellet: 


Referring  back  to  Eq.  (20)  and  Fig.  3,  the  reactant  flux  on  the 
agglomerate  core  surface  can  be  written  in  terms  of  reactant  con¬ 
centrations  using  Fields  law  for  the  ID  spherically-symmetric 
domain  as: 


Nagg  =  jj  _agg 
°2  r, 


6  /C0  -  G 


agg 


(31) 


where  Cs  can  be  expressed  in  terms  of  Co  using  Eqs.  (28)  and  (29)  as: 
Co 


Cs  = 


1 


Deff 
"  ~D~ 


(^)(Thcoth(Th)-1) 


(32) 


Similar  to  the  overpotential,  Co  is  assumed  constant  for  a  given 
agglomerate,  whereas  it  varies  across  agglomerates  within  the 
macroscale  CL  domain.  Within  the  macroscale  domain,  Co  can  be 
found  by  Henry’s  law  as: 


r  _  1  o2 
C°-~H 


(33) 


Using  the  solution  procedure  outlined  in  Ref.  49]  the  O2  con¬ 
centration  within  the  agglomerate  core  and  the  ionomer  film  (so¬ 
lutions  for  Eqs.  (22)  and  (23))  are  given  by: 


where  P0l  is  the  partial  pressure  of  O2,  and  H  is  Henry’s  constant. 
Finally,  using  Faraday’s  law  and  Eq.  (20),  the  volumetric  current 
density  is  obtained  as  follows: 


sinhfTh/^j 

C(r)  =  Cs  r  s|n^^h)  in  the  agglomerate  core 
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C0  in  the  ionomer  film 
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Table  3 

Parameters  used  in  PEM  fuel  cell  cathode  model  simulations. 


Pressure  P 

15,1987.5  Pa 

Temperature  T 

353  K 

Relative  humidity  RH 

0.5 

Stoichiometry  Ast 

2 

O2  mass  fraction  y0l 

0.209 

H20  mass  fraction  yH2o 

0.102 

GDL  porosity  fGdl 

0.4 

CL  porosity  fCl 

0.1 

Ionomer  vol.  fraction  in  agglomerate  £agg 

0.5 

GDL  permeability  fcGDL 

1.45e-ll  m2 

CL  permeability  /<CL 

1.45e-ll  m2 

Platinum  loading  mPt 

4e-3  kg  nrr2 

Platinum  particle  radius  rPt 

5e-9  m 

Agglomerate  radius  ragg 

le-6  m 

Effective  agglomerate  area  aagg 

3.6e5  m-2  m~3 

Ionomer  film  thickness  5 

80e-9  m 

Reference  02  concentration  C02jref 

0.85  mol  m  3 

Henry’s  constant  H 

31,663  Pa  m3  moD1 

Reference  exchange  current  density  io,ref 

3.86e-4  A  m~2  for  Vc  >  0.8  V 
1.46e-2  A  m-2  for  Vc  <  0.8  V 

Cathodic  transfer  coefficient  ac 

1  for  Vc  >  0.8  V 

0.617  for  Vc  <  0.8  V 

Electron  conductivity  as 

100  S  ITT1 

VI  =  — 4FV-No2  =  4FaaggN*f 


(34) 


Vi  —  4Faaggff 
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(35) 


Equation  (35)  includes  two  macroscale  variables,  Po2  and  cp 
(embedded  in  the  Th  term),  which  are  related  to  the  microscale 
parameters  C0  and  <?agg.  All  the  remaining  parameters  are  constants 
relating  to  the  material  and  geometric  properties  of  the  micro¬ 
structure.  So  Eq.  (35)  gives  the  current  generation  per  unit  volume 
of  the  CL  in  terms  of  reactant  concentration  and  overpotential  by 
representing  the  effect  of  CL  microscale  agglomerate  structure 
(without  solving  the  microscale  problem). 


2.2.  Modified  agglomerate  model 

The  modified  agglomerate  model  presented  in  Ref.  [33]  employs 
random  distributions  of  Pt  particles  inside  the  agglomerate  core 
instead  of  the  homogenous  mixture  assumption  in  the  classical 
approach.  In  the  modified  approach,  the  reaction-diffusion  problem 
needs  to  be  solved  numerically,  whereas  the  classical  model  has  an 
analytical  solution.  Therefore,  the  micro  to  macroscale  coupling 
procedure  here  is  based  on  a  numerical  solution.  In  addition,  unlike 
the  classical  spherical  ID  approach,  the  modified  approach 
employed  a  2D  cylindrical  agglomerate  configuration  for  compu¬ 
tational  efficiency  with  very  little  loss  in  accuracy.  One  can  refer  to 
[33]  for  further  details  about  the  modified  approach. 

Equation  (35)  can  now  be  applied  to  accept  numerical  results 
from  the  modified  microscale  model.  The  goal  is  to  couple  the 
modified  model  in  a  way  that  eliminates  the  need  to  numerically 
solve  for  the  microscale  model  at  every  grid  point.  Equation  (34) 
indicates  that  we  need  to  solve  for  N^g  at  the  microscale  in  order 
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to  obtain  the  local  value  of  the  current  generation  per  unit  volume 
of  the  CL.  Furthermore,  Eq.  (35)  indicates  that  N^g  can  be  expressed 
as  a  function  of  microscale  parameters  Co  and  <pagg,  which  can  be 
written  in  terms  of  macroscale  parameters  P0l  and  <p.  So  we 
simulated  the  modified  microscale  model  for  various  values  of  Co 
and  found  that  N^gg  varies  linearly  with  C0,  and  also  with  the 
macroscale  parameter  P0r  We  conducted  the  simulations  for 
agglomerate  overpotentials  <pagg  between  0  and  1  V  and  obtained  a 
curve  (called  the  agglomerate  characteristic  curve)  for  N^gg  vs.  <?agg. 
Next,  we  fitted  a  polynomial  to  this  characteristic  curve  and 
expressed  it  in  terms  of  Co.  Finally,  the  volumetric  current  density  of 
the  CL  was  obtained  as: 

Vi  =  4Faagg— (<p)  (36) 

rtL0,ref 

where  Cgg|f  =  8  mol  m-3  is  the  reference  Co  value,  and  f[cp)  is  the 
polynomial  fit  to  the  Nggg  vs.  </>agg  characteristic  curve. 

It  is  important  to  note  that  the  characteristic  curves  for  the 
modified  agglomerate  model  vary  in  the  diffusion-limited  region, 
as  shown  in  Fig.  5,  depending  on  the  specific  distribution  of 
catalyst  particles  in  the  agglomerate  core.  This  is  because  varia¬ 
tions  in  the  randomly-generated  distributions  result  in  differences 
in  diffusion  at  the  particle  level,  leading  to  different  limiting  cur¬ 
rents  as  explained  in  Ref.  [33].  Reference  [33]  provides  details  on 
the  effect  of  particle  locations  on  diffusion  losses.  Hence  it  is 
necessary  to  average  the  results  across  an  adequate  number  of 
samples  with  random  distributions.  We  investigated  samples 
containing  10,  20  and  30  random  distributions  of  catalyst  particles 
and  obtained  standard  deviations  of  their  limiting  current  den¬ 
sities  as  0.048,  0.044,  and  0.0441  A  cm-2,  respectively.  It  was 
decided  that  a  sample  size  of  20  was  adequate  to  obtain  a  reliable 
average  of  the  diffusion-loss  behavior  of  agglomerates  as  shown  in 
Fig.  5. 

We  assessed  the  accuracy  of  the  modified  agglomerate  model 
as  shown  in  Fig.  6.  First,  we  validated  our  COMSOL  code  by 
comparing  our  numerically  obtained  polarization  curve  for  the  1 D 
spherically-symmetric  agglomerate  model  with  its  analytical  so¬ 
lution.  As  Fig.  6  indicates,  the  comparison  is  perfect.  Next,  we 


Fig.  4.  2D  cylindrical  computational  domain  for  the  modified  agglomerate  model  in 
which  Pt  particles  can  be  randomly  or  selectively  distributed. 
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Fig.  5.  Agglomerate  polarization  curves  from  the  modified  approach  for  20  different 
random  particle  distributions  (mPt  =  0.4  mg  cm-2). 

investigated  the  loss  in  accuracy  incurred  by  using  a  2D  cylindrical 
agglomerate  approximation  in  place  of  a  spherically  symmetric 
one.  As  previously  mentioned,  the  2D  formulation  is  preferred  as  it 
is  computationally  more  efficient  [33].  Fig.  6  shows  that  the  po¬ 
larization  curve  from  the  2D  simulation  differs  from  the  analytical 
solution  in  the  diffusion  loss  region,  but  the  difference  is  accept¬ 
ably  small.  Next,  we  plot  the  polarization  curve  obtained  with  the 
discrete  particle  model.  It  is  immediately  apparent  that  our 
modified  particle  approach  predicts  much  higher  diffusion  losses 
compared  to  the  classical  approach.  As  explained  in  Ref.  [33],  this 
is  because  the  outermost  particles  in  a  random  distribution  may 
be  situated  some  distance  away  from  the  ionomer  film  (see  Fig.  4) 
leading  to  a  higher  effective  film  thickness,  and  hence  higher 
diffusion  losses.  Accordingly,  when  we  reduced  the  thickness  of 
the  ionomer  film  in  the  discrete  particle  model,  we  found  that  it 
was  possible  to  match  the  limiting  current  predicted  by  the  clas¬ 
sical  model.  Such  a  tuning  process  is  employed  by  many  other 
modeling  studies  which  adjust  agglomerate  structural  parameters 
to  match  experimental  results. 

2.3.  Macro-homogenous  model 

In  order  to  show  the  importance  of  multiscale  modeling,  we  will 
also  compare  the  agglomerate  model  results  with  the  macro- 
homogenous  model.  The  macro-homogenous  model  does  not 
consider  a  microscale  structure  for  the  CL;  instead,  the  CL  is  rep¬ 
resented  using  an  effective  medium  theory  for  which  it  is  only 
necessary  to  calculate  the  effective  diffusion  properties  for  charge 
and  species  transport.  The  volumetric  current  density  is  then 


Agglomerate  Current  Density  [A  cm'2] 
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Fig.  6.  Agglomerate  polarization  curves  for  various  microscale  models. 
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Fig.  7.  Mass  fraction  distribution  in  the  cathode  at  Vc  =  0.2  V  (a)  02,  and  (b)  H20.  The  cross-section  includes  the  channel,  the  GDL,  and  the  catalyst  layer. 


calculated  from  the  Butler-Volmer  equation.  The  area  term  ascaie  is 
used  as  a  fitting  parameter  as  in  Ref.  [2]  to  match  results  with  the 
agglomerate  models. 


•1  —  Of  scaled  ' 
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2.4.  Boundary  conditions  and  model  parameters 

The  model  parameters  and  operating  conditions  for  our  simu¬ 
lations  are  provided  in  Table  3.  For  the  conservation  of  mass  and 
momentum,  a  laminar  inflow  boundary  condition  was  used  with  a 
specified  inlet  velocity,  and  an  outlet  pressure  was  specified  as  the 
outflow  boundary  condition.  The  air  inlet  velocity  is  calculated 
from  the  operating  conditions  listed  in  Table  3.  The  mass  fractions 
at  the  air  inlet  are  also  calculated  from  the  given  operating  condi¬ 
tions  and  used  to  specify  the  inlet  boundary  conditions  for  the 
Maxwell-Stefan  equations.  For  the  diffusion  of  electrons,  we  set 
<p  =  0  along  the  surface  of  the  current  collector  rib,  and  zero  flux  at 
the  membrane/CL  interface.  For  ionic  charge  transport,  the  ion 
potential  at  the  current  collector  is  set  to  zero,  whereas  the  ion 
potential  at  the  membrane/CL  interface  is  set  to  ?7nco.  Symmetry 
boundary  conditions  are  applied  as  appropriate  in  the  transverse 
(x— z)  plane  (see  Fig.  2)  for  each  governing  equation. 

3.  Results  and  discussion 

The  modified  microscale  agglomerate  model  in  Ref.  [33]  was 
coupled  to  the  macroscopic  3D  cathode  model  using  COMSOL.  First, 
we  focus  on  the  general  trends  obtained  with  the  3D  single  channel 


cathode  model  using  the  modified  agglomerate  approach  (with 
adjusted  film  thickness).  The  goal  is  to  confirm  that  the  coupling 
works  properly.  Subsequently,  we  will  compare  our  3D  cathode 
model  results  coupled  using  the  modified  agglomerate  model 
approach  with  the  3D  cathode  model  coupled  to  the  classical 
agglomerate  model  and  the  macro-homogenous  model. 

Mass  fraction  contours  of  O2  and  H20  (vapor)  are  presented  for 
the  diffusion-loss  regime  in  Fig.  7a  and  b,  respectively.  As  expected, 
O2  is  being  depleted  as  it  passes  through  the  channel.  In  addition, 
the  O2  concentration  is  always  higher  under  the  channel  than  un¬ 
der  the  land,  which  implies  a  higher  reaction  rate  under  the 
channel.  Again  as  expected,  the  product  water  vapor  concentration 
increases  towards  the  outlet  of  the  channel. 

Fig.  8  shows  the  overpotential  distribution  through  the  CL 
thickness  under  three  different  operating  regimes.  Consistent  with 
the  results  in  Ref.  [19],  the  overpotential  distributions  in  Fig.  8  are 
seen  to  be  nonuniform.  As  in  Ref.  [19],  the  overpotential  is  always 
higher  under  the  land  because  of  shorter  electron  paths  that  reduce 
ohmic  losses.  The  difference  in  overpotential  between  the  land  and 
channel  regions  increases  with  reaction  rate  because  the  ohmic 
losses  become  more  noticeable  at  higher  current  densities.  It  is  also 
important  to  note  that  reaction  rate,  and  so  the  current  density 
distribution,  is  an  exponential  function  of  overpotential.  Therefore, 
even  the  small  variations  shown  in  Fig.  8  can  significantly  affect  the 
current  density  distribution. 

Volumetric  current  density  is  plotted  along  the  longitudinal  (x— 
y)  plane  occupied  by  the  membrane/CL  interface  in  Fig.  9  for  three 
different  operating  regimes.  At  low  current  densities  (NCO  =  0.3  V), 
the  reaction  rate  is  slow  enough  that  02  diffusion  is  not  the  limiting 
factor.  Therefore,  the  reaction  rate  is  maximized  under  the  land  due 
to  reduced  ohmic  losses  as  shown  in  Fig.  8.  On  the  other  hand,  at 
high  current  densities  (NCO  =  1.0  V),  mass  transport  losses  begin  to 


NCO  =0.3  [V] 


Fig.  8.  Cathode  overpotential  ( cp )  distribution  across  the  CL  thickness  from  the  modified  agglomerate  approach  along  the  transverse  (x— z)  plane  aty  =  0.95Lch.  The  bottom  boundary 
of  the  plot  corresponds  to  the  CL/membrane  interface. 
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Land  Channel  T  70576xl0‘  Land  Channel  T  1-2576xl°8  Land  Channel  ▼26837x10s 
Fig.  9.  Volumetric  current  density  along  the  membrane/CL  interface  (x-y  plane)  (A  m  3). 


dominate  and  hence  the  highest  reaction  rates  occur  under  the 
channel  region.  Finally,  at  intermediate  current  densities 
(NCO  =  0.6  V),  Fig.  9  indicates  that  the  reaction  zone  straddles  the 
junction  between  the  land  and  channel  owing  to  roughly  equal 
contributions  from  ohmic  and  mass  transport  losses;  in  fact,  the 
reaction  rate  is  maximized  at  the  corner  of  the  land.  In  addition, 
Fig.  9  also  shows  the  effect  of  mass  transport  loss  along  the  channel. 
In  all  cases,  the  reaction  rate  decreases  away  from  the  inlet,  and  the 
longitudinal  gradient  increases  with  rjuco- 

Results  presented  thus  far  show  that  the  coupled  modified 
agglomerate  approach  is  successful  at  providing  the  expected  re¬ 
sults  for  the  3D  PEM  fuel  cell  cathode.  Fig.  10  presents  polarization 
curves  predicted  by  the  classical  spherical  agglomerate  model,  the 
particle  model  with  and  without  matched  ionomer  film  thickness, 
and  the  macro-homogenous  model.  It  is  seen  that  the  macro- 
homogenous  model  cannot  produce  realistic  predictions  in  the 
mass  transport  limited  region  due  to  an  absence  of  microscale 
structure  that  could  account  for  the  dissolution  and  diffusion  of 
reactant  in  the  ionomer.  On  the  other  hand,  the  classical  agglom¬ 
erate  model  is  able  to  capture  the  sharp  drop-off  in  current  density 
in  the  mass  transport  limited  region  as  seen  in  Fig.  10.  It  should  be 
noted  that  the  ionomer  film  is  the  only  diffusion-loss  mechanism 
for  the  classical  agglomerate  which  neglects  the  significance  of 
diffusion  inside  the  agglomerate  core  [33].  On  the  other  hand,  the 
modified  particle  approach  correctly  accounts  for  the  significance 
of  diffusion  inside  the  agglomerate  core  by  considering  a  discrete 
distribution  of  particles.  Therefore,  the  particle  approach  predicts 
an  earlier  and  sharper  drop-off  of  voltage  in  the  mass  transport 
limited  region  than  the  classical  approach.  However,  Fig.  10  shows 
that  when  the  ionomer  film  thickness  is  adjusted  in  the  modified 


Current  Density  [A  cnr2] 


Fig.  10.  Comparison  of  polarization  curves  for  the  3D  cathode  obtained  with  different 
CL  modeling  approaches. 


model  to  match  that  of  the  classical  model,  the  overall  perfor¬ 
mances  are  comparable.  Hence,  for  the  comparison  purposes,  it  is 
appropriate  to  match  the  film  thickness. 

Next  we  compare  the  predictions  of  current  density  distribution 
by  the  three  models  along  a  transverse  line  on  the  membrane/CL 
interface  close  to  the  outlet.  We  investigated  the  ohmic  and  mass 
transport  loss  regions  in  Figs.  11  and  12,  respectively.  Fig.  11  in¬ 
dicates  that  for  the  ohmic  region,  all  models  predict  the  highest 
current  density  close  to  the  corner  of  the  land.  It  is  also  seen  that 
the  classical  and  the  modified  agglomerate  model  predictions  are 
close  to  each  other  throughout.  However,  while  the  agglomerate 
models  show  a  higher  current  density  under  the  channel,  the 
macro-homogenous  model  predicts  a  higher  current  density  under 
the  land.  For  the  mass  transport  region  in  Fig.  12,  the  modified 
agglomerate  model  behaves  identically  with  the  classical  approach. 
While  the  agglomerate  model  predicts  the  highest  current  densities 
under  the  channel  due  to  mass  transport  limitations,  the  macro- 
homogenous  approach  continues  to  predict  the  highest  current 
density  close  to  the  corner.  As  mentioned  earlier,  the  predictions  of 
the  macro-homogenous  approach  in  Figs.  11  and  12  are  not  realistic 
because  it  does  not  account  for  the  dissolution  and  diffusion  of 
reactant  in  the  ionomer. 

Finally,  we  investigated  the  effects  of  Pt  loading  on  the  polari¬ 
zation  curve.  Figs.  13  and  14  show  the  effect  of  Pt  loading  on  po¬ 
larization  curves  as  predicted  by  the  classical  and  modified 
agglomerate  models,  respectively.  Although  the  trends  are  similar 
for  both  models  in  the  activation  and  ohmic  loss  regions,  the  pre¬ 
dictions  diverge  in  the  mass  transport  region.  The  classical 
agglomerate  model  predicts  nearly  identical  limiting  currents  as  in 
Ref.  [19]  despite  a  fourfold  increase  in  Pt  loading  as  seen  in  Fig.  13. 


Position  (land  to  channel  at  z=0)  [mm] 

Fig.  11.  Comparison  of  current  density  distribution  from  different  CL  modeling  ap¬ 
proaches  along  a  line  at  the  membrane/CL  interface  at  y  =  0.98LCh  and 
Lvg  =  0.4  A  cm-2. 
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Position  (land  to  channel  at  z=0)  [mm] 

Fig.  12.  Comparison  of  current  density  distribution  from  different  CL  modeling  ap¬ 
proaches  along  a  line  at  the  membrane/CL  interface  at  y  =  0.98LCh  and 
4vg  =  0.77  A  cm-2. 


Fig.  14.  Polarization  curves  for  various  Pt  loadings  (mPt)  in  mg  cm  2  predicted  by  the 
modified  agglomerate  model. 


This  is  because  diffusion  losses  in  the  classical  agglomerate  model 
are  dominated  by  the  ionomer  film  which  imposes  the  same 
limiting  current  irrespective  of  Pt  loading  inside  the  core  [33].  On 
the  other  hand,  experimental  results  in  Ref.  [50]  show  that  Pt 
loading  affects  the  performance  and  the  limiting  current.  The 
modified  approach  correctly  accounts  for  the  significance  of  diffu¬ 
sion  inside  the  agglomerate  core  by  modeling  the  particle-level 
effects.  Therefore,  the  particle  model  is  sensitive  to  Pt  loading  in 
the  mass  transport  loss  region  as  shown  in  Fig.  14. 

It  should  be  noted  that  it  is  possible  obtain  different  limiting 
currents  even  with  the  classical  agglomerate  model  by  setting  re¬ 
lations  between  the  agglomerate  and  CL  structural  parameters.  For 
instance,  Refs.  [20-22,26,27]  incorporated  changes  in  the  CL 
structure  by  varying  the  number  of  agglomerates  per  unit  volume 
which  concurrently  changes  the  total  effective  agglomerate  surface 
area  (aagg).  However,  none  of  these  is  equivalent  to  incorporating 
the  effect  of  particle-level  diffusion  losses  in  the  microscale  model. 

4.  Summary  and  conclusions 

.In  this  study,  we  employed  a  modified  agglomerate  approach 
based  on  discrete  catalyst  particles  [33]  to  build  a  3D  model  for  the 
PEM  fuel  cell  cathode.  The  effect  of  random  agglomerate  structure 
and  nanoscale  diffusion  losses  are  incorporated  within  a  3D  PEM  fuel 
cell  model  for  the  first  time.  The  microscale-to-macroscale  coupling 
procedure  employed  in  the  classical  agglomerate  approach  in  Ref. 
[19]  was  modified  to  permit  the  use  of  numerical  results  from  the 


Fig.  13.  Polarization  curves  for  various  Pt  loadings  (mPt)  in  mg  cm  2  predicted  by  the 
classical  agglomerate  model. 


modified  agglomerate  microscale  model.  The  novel  3D  cathode 
model  was  found  to  produce  expected  results  for  species  concen¬ 
tration,  overpotential,  and  reaction  rate  distributions  in  the  cathode. 

By  comparing  the  predictions  of  the  modified  model  with  the 
classical  and  macro-homogeneous  models,  it  was  shown  that  the 
modified  approach  provided  results  that  matched  with  the  classical 
approach  under  all  operating  regimes.  However,  the  macro- 
homogenous  model,  which  ignores  the  importance  of  microscale 
effects,  fails  to  successfully  predict  mass  transport  losses. 

In  addition,  the  modified  approach  is  superior  to  the  classical 
approach  in  predicting  the  effect  of  Pt  loading  inside  the  agglom¬ 
erate  core.  The  modified  model  exhibits  sensitivity  to  Pt  loading 
even  in  the  mass  transport  loss  region,  whereas  the  classical 
approach  provides  nearly  identical  results  irrespective  of  Pt 
loading. 
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Nomenclature 

Ao  surface  area  per  unit  weight  of  Pt  [cm2  g-1  ] 
aPt  total  effective  reaction  area  [m2  m-3] 
aagg  effective  agglomerate  surface  area  [m2  m-3] 

C02  local  oxygen  concentration  [mol  m-3] 

C  oxygen  concentration  in  agglomerate  core  [mol  m-3] 

Cf  oxygen  concentration  in  ionomer  film  [mol  m-3] 

Cs  agglomerate  surface  oxygen  concentration  [mol  m-3] 

Co  ionomer  film  surface  oxygen  concentration  [mol  m-3] 
Q)2,ref  reference  oxygen  concentration  [mol  m-3] 

Dij  binary  diffusivities  of  i  and  j  [m2  s-1] 

Dfjf  effective  binary  diffusivities  of  i  and  j  [m2  s-1] 

D  diffusivity  of  oxygen  in  ionomer  [m2  s-1] 

Deff  effective  diffusivity  of  oxygen  in  ionomer  [m2  s-1] 

Eth  Nernst  voltage 

H  Henry’s  constant  [atm  m3  mol-1] 

i  current  density  [A  cm-2] 

V.i  volumetric  current  density  [A  cm-3] 
if  reference  exchange  current  density  [A  cm-2] 
kc  reaction  rate  constant  [s-1] 

Mn  molecular  weight  of  gas  mixture  [kg  mol-1] 
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Mj  molecular  weight  of  species  j  [kg  mol-1  ] 
mpt  Pt  loading  [mg  cm-2] 

Nq2  oxygen  flux  [mol  m-2  s-1] 

P  pressure  [atm] 

Pq2  oxygen  partial  pressure  [atm] 

Ph2  hydrogen  partial  pressure  [atm] 

Pt|C  platinum  carbon  weight  ratio 

R0l  volumetric  oxygen  reaction  rate  [mol  m-3  s-1] 

PHo2  volumetric  reaction  rate  for  water  vapor  [mol  m-3  s-1] 
ragg  agglomerate  radius  [m] 

rPt  Pt  particle  radius  [m] 

T  temperature  [I<] 

Th  Thiele  modulus 

tci  catalyst  layer  thickness  [m] 

ac  charge  transfer  coefficient 

«od  osmotic  drag  coefficient 

d  ionomer  film  thickness  [m] 

£gdl  porosity  of  GDL 

£Cl  porosity  of  CL 

£agg  ionomer  volume  fraction  in  agglomerate 

Vnco  nominal  cathode  overpotential 

A  local  water  content 

It  dynamic  viscosity  [kg  m-1  s-1] 

p  density  [kg  rrT3] 

as  electron  conductivity  [S  m-1] 

offf  effective  electron  conductivity  [S  m-1] 

am  ionic  conductivity  [S  m-1] 

effective  ionic  conductivity  [S  m-1] 

<p  local  overpotential  [V] 

0 s  electron  potential  [V] 

0m  ionic  potential  [V] 

0rh  local  relative  humidity 

0  Thiele  modulus 

oji  mass  fraction  of  species  i 
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